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INTRODUCTION 

The  geometry  for  the  following  analysis  is  shown  in  Figure  1.  A  point  source  of  current  / 

is  located  on  the  z  axis  at  z  coordinate  h.  Space  is  divided  into  three  regions  with  electrical 
conductivities  as  follows: 


o  =  0, -oo < z  ^0 

(region  0) 

da) 

a  =  c1,0<z£d 

(region  1) 

(1  b) 

<5  =  <5 2,d  <z<°° 

(region  2) 

(1  c) 

region  0  _ 

surface 


FIGURE  1.  GEOMETRY  FOR  POINT  CURRENT  SOURCE  FIELDS 
The  source  is  located  in  the  middle  region.  It  has  a  potential  of  the  form 


<*>,= 


47tC7i ^jr2+(z  -h)2’ 


which  has  the  integral  representation 


oo 

®>-4 


Elsewhere,  the  potential  function  satisfies  Laplace’s  equation 

V2O(r,e,z)  =  0, 


(2) 


(3) 


(4) 
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which  has  the  axially  symmetric  solution  in  cylindrical  coordinates 

oo 

Ow=  J  d\J0{Xr)[a{X)eu  +  b{X)e-^l 

o 

The  potential  in  the  three  regions  then  has  the  physically  acceptable  forms 


and 


^>0  =  J dXJ0(kr)a0(k)eXi  (-«  <  z  < 


0) 


oo 

*1  =  J  dXJ0(kr)\ ^e-^-v+a^  +  bWe-* 

o  L  1 


oo 

<&2  =  J  J0(^r)b2(K)e~^  ( d<z<  °°). 


(0  <  z  <  c?) 


BOUNDARY  CONDITIONS 

The  boundary  conditions  require  continuity  of  the  normal  component  of  the  current 
and  the  tangential  component  of  the  electric  field  at  each  boundary.  The  electric  field  and 
density  are  given  by 

E  =  -V<b  and  J  =  gE. 

Continuity  of  the  normal  current  density  at  each  boundary  leads  to  the  equations 

a1(X)-bl(k)  =  --^-e^ 

4710! 


and 


a  A*”  ~  bA)e~U  +  ebA^  =J—e~Hd-h) 

4k  Gj 


where 


e  =  G2/Gj. 

Continuity  of  the  tangential  electric  field  at  the  bottom  gives  the  equation 
aA)eU  +  bx(K)e~U  -  b2(k)e-u  = 

47CGi 


(5) 

(6) 

(7) 

(8) 

density 

current 

(9) 

(10) 

(11) 

(12) 

(13) 
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Equations  (11)  ,(12),  and  (13)  are  sufficient  to  determine  a,(X,) ,  fc,(X) ,  and  b2(k).  This  allows 
the  computation  of  all  fields  in  the  middle  and  lower  regions,  which  are  of  primary  interest  here. 
The  function  a0(X)  can  be  determined  from  continuity  of  the  tangential  electric  field  at  the  air- 
conductor  interface,  if  needed. 


With  the  introduction  of  the  parameter  Q  =  (1  -e)/(l  +e),  the  solutions  for  the  coefficients 
in  the  upper  and  lower  conducting  regions  can  be  written 


and 


-nj 


a1(k)  = 


471©!  1-Qe 


Qe'^n(e»  +  e^) 


bA)  = 


4ko1 


r.  -ild 


1-Qe 


(14) 

(15) 


b2(X)  = 


I  (1-0)0?**  +  *-**) 
47to2  1  -Qe~^ 


(16) 


The  parameter  Q  has  the  range  of  values  -1  <  Q  <  1,  where  0  =  1  corresponds  to  region  2  being  a 
perfect  insulator,  0=0  corresponds  to  regions  1  and  2  being  identical,  and  Q  =  -1  corresponds  to 
region  2  being  a  perfect  conductor. 


POTENTIAL  AND  CURRENT  DENSITY 


The  potential  for  the  middle  and  lower  regions  is  given  by 


0l  =4w  J dXjo(^){e~Mz~h'  +e~Ht+h)  (17) 

1  o 


— [cw.«)+ e-»<=**)+ ev.-*)+ ,, 

l  ^ 


and 


eo 


V  e-Kz-h)+e-Kz+hy 


1-Qe 


-2 Xd 


(18) 


Integrals  involving  the  bessel  functions  are  notoriously  difficult  to  perform  numerically. 
However,  a  geometric  series  expansion  valid  for  |  Q\  <1 


1 

1-Qe-** 


£  0  V2'’*'' 


(19) 
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and  the  integral  representations  exemplified  by  Equations  (2)  and  (3)  can  be  used  to  write  the 
potential  in  the  upper  and  lower  conductors  as 


1 


l 


+  I  Qn[- 


1 


and 


4nai  Vr2  +  (z -hf  yjr2+(z  +  h)2  "  =  >  ^r2  +  (2nd -z -hf 

1  _ 1 _ ^ _ 1 _ ni 

> lr2  +  (2nd  +  z  +  h )2  yjr2  +  (2nd -z  +h)2  ylr2  +  (2nd+z~hf ^ 


<^=— (l-G){  1  1 


(20) 


=- — (l  —  8){~r  +-r 

47102  yr2  +  (z  -hf  yjr2  +  (z+h)2 


(21) 


+  1  Q‘ 


n  =  1  l^lr2+(2nd  +  z-h)2  ^r2  +  (2nd  +  z  +h)2 


}• 


By  combining  Equation  (9)  with  Equations  (20)  and  (21),  expressions  for  the  cartesian  components 
of  the  current  density  in  the  middle  and  lower  regions  can  be  derived.  They  are 


Jlx(y)~ 


x(y)I 


|  ,  ^  1 
471  > l(r2  +  (z-hff  yl(r2  +  (z+h)2f  n  =  1  V (r 2  +  (2nd  -z-h ff 

1,1.1 

V (r2  +  (2nd  +  z+  h )2)3  'l (r2  +  (2nd  -  z  +  hff  + ^ (r2  +  (2nd  +  z  -  hff^ 


(22) 


Jiz=-7Z{ 


z+h  +  -  ^  2**-z-A 


4rc  *N l(r2  +  (z-hff  ^l(r2+(z+hff  "  =  1  yj (r2  +  (2nd  -  z  -  hff 


(23) 


2nd  +  z+h 


2nd  -z+h  _  2nd  +  z-h 

+  , - - - ]}. 


V  (r2  +  (2nd  +  z+  h  )2)3  V(r2  +  (2  nd-z+hff  ^l(r2  +  (2nd  +z -hff 
x(y)1  ..  _  -  _J  1  1 


- 


47C 


-(1-0)  I  0" 


"=°  LV(r2  +  (2nJ  +  z-/i)2)J  '\j(r2  +  (2nd  +  z  +  hff  _ 


(24) 


and 


47C  n  =0 


2/t<f  +z-h 


2nd  +  z+h 


_^(r2  +  (2nd  +  z-hff  (r2  +  (2nd  +  z  +  hff 


(25). 
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MAGNETIC  FIELD  DUE  TO 
DISTRIBUTED  CURRENT 

The  magnetic  field  of  the  distributed  current  can  be  constructed  by  means  of  a  careful 
application  of  Ampere’s  Law,  applied  to  a  properly  constructed  system;  that  is,  one  in  which  charge 
conservation  holds.  To  accomplish  this,  the  point  current  source  is  imagined  to  be  fed  by  an  insulated 
wire  that  enters  region  1  perpendicular  to  die  boundaries,  and  extends  into  medium  0  for  a  great 
distance,  ultimately  returning  to  a  point  current  sink  in  region  1  at  a  great  lateral  distance  from  the 
current  source.  Application  of  the  Biot-Savart  Law  to  the  wire  shows  that  if  the  wire  extends  far 
enough  from  the  boundary  in  the  z  direction  and  returns  to  the  medium  sufficiently  far  away,  the 
magnetic  field  of  the  wire  segment  is  indistinguishable  from  that  of  a  semi-infinite  wire  extending 
in  the  -z  direction.  In  addition,  if  the  current  sink  is  sufficiently  far  away,  its  fields  are  negligible 
at  the  current  source,  and  axial  symmetry  about  the  incoming  wire  can  be  applied.  This  is  illustrated 
in  Figure  2. 


*  FIGURE  2.  MAGNETIC  FIELD  OF  WIRE  AND  DISTRIBUTED  CURRENTS* 


Ampere’s  Law  can  now  be  applied  to  the  current  distribution.  Ampere’s  Law  takes  the  form 

jn-d\  =  jj'dA,  (26) 

which  states  that  the  line  integral  of  the  magnetic  field  around  a  closed  path  is  equal  to  the  integral 
of  the  normal  current  density  over  the  enclosed  surface;  that  is,  the  net  current  through  the  loop.  In 
the  present  application,  symmetry  about  the  wire  axis  can  be  used  to  write 

j?L>dl  =  2KrH9  (27) 


“Note  that  the  magnetic  field  is  entirely  horizontal,  there  is  no  vertical  field  from  the  wire  current  or  the  distributed 
currents  from  the  point  source. 
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where  r  is  the  radius  of  a  circle  centered  on  the  wire  axis  and  /7e  is  the  tangential  component  of 

magnetic  field  on  the  circle,  which  is  the  same  everywhere,  independent  of  0.  Now,  if  the  integral 
on  the  right  hand  side  of  Equation  (26)  can  be  evaluated,  the  magnetic  field  in  regions  1  and  2  can 
be  determined  directly .  To  accomplish  this,  region  2 ,  d  <  z  <  is  first  examined,  with  reference 
to  the  notation  of  Figure  1.  The  right-hand  rule  is  used  to  assign  the  sense  of  the  currents  and  fields, 
and  then 


dA 


-J 


J2idA  = 


2  k  r  r 

j  dQ  J duuJu  =  2tc  J duuJ2z 


(28) 


where  J2z  is  independent  of  0. 

Substitution  of  Equation  (25)  into  Equation  (28)  and  application  of  Equations  (26)  and  (27) 

yields 


^2e  =  4^{2~(1~g)(  /  2Z  h  2+l  >Z+h 

4nr  Vr2  +  (z -h)2  V r2  +  (z+hf 


+  lQn\ 

n  =  1 


2nd  +z-h  2nd  +z+h 

I  -H - 7"  ' 

^r2  +  (2nd  +z -hf  yjr2  +  (2nd  +  z  +h)2 


)}• 


(d  <z  <° °) 


(29) 


Repeating  this  analysis  for  region  1,  for  h  <  z  £  d ,  and 
(25), 


using  Equation  (23)  instead  of  Equation 


'  -  1  [o  z~h  z+h  |  -  2nd  —z—h 

4nr  y/r2  +  (z-h)2  y/r2  +  (z+h)2  »  =  >  ^jr2  +  (2nd -z  -hf 

_ 2nd  +z+h  |  2nd  -z+h  2nd  +z  -h  11 

ylr2  +  (2nd+z+h)2  ylr2  +  (2nd-z +h)2  ^  r2  +  (2nd  +  z  ^hf]} 

(h  <z< d). 


(30) 


Next,  region  1  is  examined  for  the  case  0  <  z  <  h .  Consistent  application  of  Ampere’s  Law 
for  this  case  requires  that  the  contribution  of  the  line  current  be  included,  and  then 


2nrHe  =  I  +  2k 


j  duuJz 
o 


(31) 


and  this  time  substitution  from  Equation  (23)  with  careful  attention  to  the  relative  values  of  z  and 
h ,  yields 
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r  _  1  [2  z~h _ z+h  |  £  2 nd-z-h 

47Cr  Vr2  +  (z  -/i)2  ^r2+(z+h)2  n=1  yjr2  +  (2nd  -z-h)2 

2nd+z+h  ^  2nd -z+h  2nd  +z  —h 

^r2  +  (2nd  +z  +  hf  yjr2  +  (2nd -z+h)2  yjr2  +  (2nd  +  z -h)2 

(0£z  <h). 


which  is  identical  in  form  to  Equation  (30),  and  the  magnetic  field  is  continuous  at  the  electrode 
(z=h)  plane,  as  it  must  be  to  be  physical. 

Now,  to  isolate  the  magnetic  field  due  to  the  electrode  distributed  currents  alone,  the  Biot- 
Savart  magnetic  field  of  the  wire  alone  is  constructed  and  it  is  subtracted  from  the  expressions  in 
Equations  (29),  (30),  and  (32).  Then 


(z-h)  ' 

V r2+(z-h)2 

and  the  electrode  magnetic  field  alone,  everywhere  in  region  1,  is  given  by 

=  —  - - +  £  . 

4nr  V r2  +  (z+h f  "=1  yjr2  +  (2nd -Z-h)2 

2nd  +  z  +  h  |  2nd -z+h  2nd  +z~h 

yjr2+(2nd  +z+h)2  yjr2+(2nd -z+h)2  ^r2  +  (2nd  +  z-h)2 

(0  <z<d), 


0  47tr 


(33) 


(34) 


and  in  region  2,  the  field  is 
f  L .  z-h 


H‘e=  A 
20  47tr 


i+-~r2-  -(i-g)Ifil 

Vr2  +  (z-/i)2  n=0 


2nd  +z  —h 


2nd  +  z  +h 


y lr2  +  2nd  +z  -h2  yjr2  +  2nd  +  z+h2 


"•  (35) 


Then  the  cartesian  components  in  either  region  are  given  by 


and  He=-Hl. 
yr 


(36) 


NON-CONDUCTING  BOTTOM  (Q=l) 


Region  1 


In  the  limit  of  a  non-conducting  bottom,  Q  =  1 ,  and  the  series  solutions  for  the  fields  diverge. 
To  deal  with  the  fields  in  region  1,  the  integral  expressions  defining  them  will  be  revisited.  The 
current  density  is  determined  from  Equations  (9)  and  (17)  and  has  the  form 
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and 


«o 

+  (37) 

e-M 

+  l-e-^ 


4  =£  J dUJ0(h-)[sSnh  -h)e-',‘-'"+e-H•'h,  (38) 

0 

-Ud 

+l~^[-eK^h)+e-*l+h)-eKt-h)+e-*-hy]}. 


[eMz *h) + e~Uz+h)  +  eMz~h)+  e~Hz _*>1} 


The  case  z  >  h  is  explicitly  selected.  The  end  result  is  the  same  for  any  case.  The  expression 
for  the  electrode  magnetic  field  is  then  obtained  by  combining  Equations  (26)  through  (28),  (33), 
and  (38).  The  result  is 


He10  =  -/~(- 1+-,  Z  h  +r 
4 ™  V r2  +  (z-h )2 


«o 


eHz-h)  +  e-Mt+h) 


(39) 


-7U 


1-e 


-2Xd 


In  each  of  Equations  (37)  through  (39)  a  common  denominator  is  used  and  terms  are  combined 
into  hyperbolic  functions.  Then,  the  fields  have  the  forms 


/  7*(y) 

■/uw~  4rcr 


f  dXU ,  (Xr)  1  -C0Sh[W  - 1  +  *  )1  +  coshim  -  z  - 
J  1  1  sinh(M) 


h)\ 


(40) 


J  =—  f  dXKJ  (\r)[  sinh^(^  ~z  +  ft)]  +  sinh[X(J-z  -h)\ 
u  4n  J  ol  ;{  sinh(XJ) 


(41) 


and 


*  h 


+r 


OO 

J  dXJ.iXr) 


V r2  +  {z-hf 

sinh[A.(t/  -z+h)]  +  sinh[A,(J  ~z~h)] 
sinh  (Kd) 


)• 


(42) 


There  are  three  integral  types;  they  are 
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and 


ii  = 


cosh(X^) 

sinh(X^) 


J  dXkJ0(kr) 

0 


sinh(Xa) 

sinh(Az?) 


^3  = 


J  dXJx(kr ) 


sinh(Azz) 
sinh(X^)  ‘ 


(43) 


(44) 


(45) 


Contour  Integration-Residue  Theorem.  If  the  integrands  of  ij,  i2,  and  t3  are  analytically 
continued  into  the  complex  plane  X  ->  w  =  u  +iv,  the  denominator  sinh(W)  will  have  zeros  for 
wd  =  inn,  where  n  =  0,±1,±2,  •••.  This  suggests  that  the  required  integrals  be  related  to  integrals 
over  a  closed  contour  containing  these  zeros,  and  the  Cauchy  residue  theorem  be  exploited.  An 
appropriate  contour  is  shown  in  Figure  3.  The  contour  is  closed  in  the  upper  half  plane,  and  the 
zeros  involved  are  those  for  n  >  0. 


The  integrals  of  interest  are  taken  over  the  interval  I2  in  Figure  3,  in  the  limit  that  la  goes  to 
0,  and  I2  becomes  infinite.  This  integral  will  be  related  to  the  closed  contour  integral  over  the 
successive  intervals  Iu  Ib,  I2,  and  Ic.  To  this  end,  the  integral  over  I2  will  be  expressed  in  terms  of 
an  integral  over  both  Ix  and  /2.  To  accomplish  this,  bessel  functions  J„  will  be  expressed  in  terms 
of  hankel  functions.  This  relationship  is  analogous  to  the  representation  of  the  cosine  function  by 
means  of  complex  exponentials,  which  is  a  trick  used  in  fourier  transform  theory  to  extend  integrals 
from  the  positive  real  axis  to  the  entire  real  axis.  The  hankel  functions  are  defined  by 
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^(zW.W  +  iT.fe)  (46) 

H?\z)=Jn(z)-iYn(z) 

where  Yn  is  the  bessel  function  irregular  at  the  origin.  The  bessel  function  Jn,  which  is  regular  at 
the  origin,  now  can  be  expressed  as 


4fe)=f[»“(z)+Hf)(z)].  (47) 

Next,  it  is  noted  that  the  hankel  functions  have  the  symmetry  properties 

H(0'X-z)  =  -H^(z)  (48) 

H[lX-z)  =  Hl2Xz). 

Inserting  Equation  (47)  in  Equation  (43),  for  example,  and  using  the  symmetry  properties  of  all  the 
factors  in  the  integrand,  gives 


(dXU^Xr) 

h 


cosh(Xa) 

sinh(A/7) 


1 

2 


cosh(Xa) 

sinh(Xd) 


cosh(Azt) 

sinh(A/7) 


1 

2 


h  -u 


cosh(-An) 

sinh(-A/f) 


/,&/ j 


cosh(Xa) 

sinh(A/7) 


(49) 


where  the  relation  has  been  used  for  the  intervals  —  (— 72)  =  7j.  Thus,  the  integral  over  72  is  replaced 

by  an  integral  over  7,  and  72,  with  Jx  replaced  by  H^/2.  A  similar  operation  can  be  done  for  the 
integrals  in  Equations  (44)  and  (45). 


Now,  using  a  shorthand  notation,  relations  between  the  integrals  over  the  various  contour 
segments  can  be  written,  where  it  is  understood  that  the  equations  hold  in  the  limit  that  the  inner 
semicircle  shrinks  to  zero,  and  the  outer  semicircle  recedes  to  infinity.  Then 


(50) 


and 
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H+M 

4*4  h  4 

where  $  is  the  integral  around  the  closed  contour.  Thus 

4  4  4 

The  following  generic  results  hold  for  the  three  integrals  of  interest. 

f— >2t o'  isn 

J  n  =  l 

J  -»  0,  (all  cases) 


(51) 


(52) 


(53) 

(54) 


J  -4  0,  (all  cases) 


(55) 


J  ->  (simple  limit) 


(56) 


Here,  the  S„  are  the  residues,  defined  by 


(dwf(w)~  ((ivt'— —————— —s  i dw 

J  J  W-W0  J 


g(w) 

—  —  Y  uw - , 

W-Wo  J  W-W0 

Sn  =  lim  g(w)  lo=inxld . 


(57) 


To  establish  the  limiting  values  of  the  various  integrals,  the  limiting  forms  of  the  functions 
involved  are  examined; 


70(w)-»l,  7j(«)->m,  m -» 0  (58) 

cosh(w)->l,  sinh(vv)  ->w,  |w| ->0  (59) 

^V)^-Dn(w/2)+y],  J?<V)->— ,  |w| ->0  (60) 

7t  KW 

(Y  =  0.5772156649--  -Euler’s  constant) 
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and 


^oV)-»— ■ ,  #i(V)  ->  |  vv|  — > °°.  (61) 

Vtch'  yiTw 

Applying  Equations  (58)  and  (59)  to  the  integrands  of  i„i2  and  t3  shows  that  the  integrands 
approach  0  as  X  approaches  0,  establishing  Equation  (54).  For  the  integrals  on  contour  Ic,w  =Re‘*, 
and  the  hankel  function  factor  in  the  integrands  has  the  behavior 


e~rR  tin* 

i —  >  0,  R  — > 

yjrR 


(62) 


The  other  factors  are  all  proportional  to  the  limiting  form 


(63) 


For  <|>  =  n/2,  the  associated  factors  are  purely  oscillatory,  and  are  dominated  by  the  decaying  hankel 
function  factor.  For  other  angles 

a-d  =d-z+h-d  =  -(z-h)<  0,  or  a -d  =d -z-h -d  =-{z+h)<  0  (64) 

and  the  factor  vanishes  as  R  ->  °o.  This  establishes  Equation  (55). 

To  determine  the  limiting  forms  for  the  integrals  over  the  contour  4,  let  w=  he  and  explicitly 
list  the  integrals.  They  are 


0 


2i>  _ 4 


n hre'*hde'* 


o 


(65) 


r  ( 
In 
L  V 


M+Y 


hae 


if 


hde 


>♦ 


0,  h  — ^  0 


(66) 


and 


1  f  if  hae'*  a  a 

.  dtyihe  —  dty — > — h  — >0. 

2d  7t hre  *hde  *  Tird  J  rd 


(67) 


The  residues  are  constructed  using  a  simple  limit  process,  as  indicated  in  Equation  (57).  Some 
useful  expressions  for  the  evaluations  are 


sinh(dw)  =cosh(dwn)d(w  -wj  +  •••,  w  wn 
where  dw„  is  a  zero  of  sinh(dw), 


(68) 
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cosh(/  u)  =  cos(n  )  (69) 

sinh(zw)  =  i  sin(w) 


and 


Ho\iu)  =  -—K0(u )  (70) 

71 

H?Xiu)  =  --Kl(u) 

7C 


where  the  Kn  are  the  modified  bessel  functions  regular  at  infinity.  Using  these  results,  the  residues 
can  be  written  as 


i-rrd), . ,COSh(flVV) 


-wH?’(rw) 


inn , 


(w  -inn/d)  -4  '-^HfXinnrtd) 


sinh(zfu') 


cosh(zn7ta/£() 

~HTd 


(71) 


=  n(-lfKx(nnrld)cos(nnald ) 


1 

2 


rwfff(rw) 


sinh(aw) 

sinh(dw) 


Z/ITt, 


(w  -  inn/d)  — » vnnr/d) 


sinhji  nna/d) 
(-l)V 


(72) 


J 

=— n(-l)  K0(nnr/d)sin(nna/d ) 


and 


1  sinh(aw)  .  1  smh(inna/d) 

-H\  \rw)  .  :  (w  -  inn/d)  -> -H\'(innr/d) - 1 - - 

2  1  sinh(<fw>)  2  1  (-1  yd 


=  — -(-If  Kx(nnr/d)  sin(nna/d). 
nd 


(73) 


Now  the  explicit  forms  for  the  fields  can  be  constructed.  In  doing  so,  the  arguments  of  the 
trigonometric  functions  are  expanded  and  integral  multiples  of  7t  are  eliminated.  Now 


i  = 


(74) 


and  Equations  (40)  through  (42),  Equation  (53),  Equations  (65)  through  (67),  and  Equations  (71) 
through  (73)  are  combined.  The  results  are 


Ix(y ) 
2ndr 1 


1+—  X  nK1(nnr/d)[cos{nn(z-h)/d}  +  co&{nn(z+h)/d}] 

y  Cl  n-1 


(75) 
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4 i=—2  £  nK0(nnr/d)[sin{nK(z-h)/d}  +  sin{nii(z+h)/d}] 

2 a  <1  =  1 


(76) 


and,  subtracting  the  Biot-Savart  field  (Equation  (33)), 

ut  I  z-h  2 z 

(77> 

Ttr  " 

+-J  X  4(«rcr/tf)[sin{n;t(z  - A)/J} +  sin{/i7t(z  +h)/d}]). 

u  n  =  1 

Region  2 

Inspection  of  Equations  (24),  (25),  and  (35)  shows  that  the  fields  reduce  to  the  simple  forms 

(78) 

(79) 


4 x(y)  ~  0 


4=0 


and 


J4=- 


47tr 


1  + 


z -h 


y lr2  +  (z -hf 


(80) 


PERFECTLY  CONDUCTING  BOTTOM  (Q=-l,  REGION  1  ONLY) 


Region  1 


In  this  case,  when  the  various  terms  in  region  1  are  combined  to  give  hyperbolic  functions, 
the  resulting  integrals  are 


oo 

<  _7x(y)  f  n,  f  sinh[X.(J -z  +/i)]  +  sinh[>,(J -z -h)] 

UM  4nr  J  '[  cosh(M) 


cosh[X(^  -  z  +  h )]  +  cosh[^,(J  -z-h)] 
cosh(ta0 


and 


cosh{7.(J  -  z  +  h )]  +  cosh[A,(J  -z-h)] 
cosh  (Kd) 


The  three  types  of  integral  now  are 


(81) 


(82) 


(83) 
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J  dXU^lr) 

o 


sinh(Xg) 

cosh(X^) 


(84) 


h  = 


J  dXkJQ(kr) 

o 


cosh(Xa) 
cosh  (Ad) 


(85) 


and 


h  = 


J  dUt(Xr) 


cosh(Aa) 
cosh(Ad)  ’ 


(86) 


Again  the  construct  of  Equation  (49)  can  be  used  to  introduce  the  hankel  function,  and  analytic 
continuation  and  the  residue  theorem  can  be  used,  where  now  the  singularities  occur  for 

z=i(n  +  l/2)n/d,  n  =0,±1,±2,...  (87). 

Again, 

i=fJ-M  <88> 

4  4  4 

and,  again 

/„->  0  and  Ic  -» 0.  (89) 

Now,  for  4 


2  2i>  -2/  a/te'* 
Khre'*  1 


»0 


(90) 


2/ 

7C 


In 


Are'* 


A  n 


+y 


(91) 


1 

2 


0 


J*  dtyihe 

K 


»  ~2  i  1  .  1 
7t/zrg'*l  r 


(92) 


and  for  the  residues 
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rH'Xrw)^kw-i(-n+my,ild] 


•  -[/ (n  + 1  l2)nld]H{'\Hn  +  H2)%rld]  sinh[i(n  +  \l2)nald] 


i( 


i(n  + 1/2)  (-1) 


n  + 1 


-^[(n  +  \l2)Krld]s\n[{n  +  \l2)nald} 


+ m)m 


•-[(n  + 1  l2)nld]H(o\i{n  +  l/2)itr/d]cosh[i(n  +  1/2  )na/d] 


(- 


i(n  +  l/2)(-l)"  +  1 

-  Ti - Ko[(n  +  l/2)nr/d]cos[(n  +  1/2)7 za/d] 


Then,  using 


lw'<1W)ifS[,v-,'(”+1/2),'W] 

1  x^i0)[i(«  +  l/2)7tr/J]cosh[/(n  +  1/2)71; aid] — - — 
2  i(-l)nd 

/(— 1)" 

=  ~Kd  XiKn  +  l/2)nrld]cos[(n  +  l/2)na/d]. 


■-h- 


2m  15.-/, 


n  =  1 


gives 


Ix(y)  00 

Ju{y)  =  T-~  l(n  +  \!2)K,[{n  +  1/2)tc rid] 

Ird  »-o 

{cos[(n  +  l/2)7t(z -h)/d]  +  cos[(/i  +  1/2)7 i(z+h)/d]} 


Ju=—r2l(n  +  l/2)K0[(n  +  1/2)7 zr/d] 

2a  «  =o 

{sin[(n  +  l/2)7t(z  -h)/d]  +  sin[(/i  +  l/2)7t(z  +h)/d]} 
and,  again  subtracting  the  Biot-Savart  field  (Equation  (33)), 


1 

-i  yd 


i 

i  yd 
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Hie  =  J~(l  +  ,  *  k  7— T  S  *iK»  +  m*r/d] 

4%r  ylr2  +  (z-h)2  d  -=° 

{sin[(n  +  1/2)tc(z  -/i)/^]  +  sin[(n  +  l/2)7t(z  +/i)/d]}). 

Region  2 

In  region  2,  the  integral  forms  for  the  fields  are 


(99) 


,  x(y)I 
J™~  2k r 


7.  f  ,,-Hz+h),  -Ki-hy 

j  ^i(^)  ~ 


(100) 


ic  r 

0  L 


-2X4 


(101) 


and 


«-4^{-1  + 


z -h 


e» 

+2r  J  JXJ^Xr) 


slr2  +  (z-h)2 

e-Hz+h)+e-Mz-hy 


1  +  e 


-2X4 


}• 


(102) 


In  these  integrals,  the  denominator  can  be  converted  to  a  hyperbolic  cosine,  but  the  numerators 
will  be  linear  combinations  of  hyperbolic  sine  and  cosine  functions.  Some  of  the  terms  can  be 
rewritten  using  a  construct  similar  to  Equation  (49),  and  can  be  treated  with  contour  integration 
techniques.  However,  the  remaining  terms  do  not  have  the  appropriate  symmetry  to  accomplish 
this.  These  remaining  terms  have  symmetry  such  that  the  integral  can  be  extended  over  the  real 
axis  while  retaining  die  bessel  function  form,  and  there  is  no  singularity  at  the  origin  as  there  is 
with  the  hankel  function.  However,  when  the  integrands  of  these  terms  are  analytically  continued 
into  the  complex  plane,  the  bessel  functions  diverge  in  either  half  plane,  and  the  outer  contour  cannot 
be  closed.  Thus,  there  is  no  way  to  explicitly  evaluate  these  integrals.  As  this  particular  case  is 
more  of  a  curiosity  without  practical  significance,  it  will  not  be  pursued  further. 


DISTRIBUTED  CURRENT  MAGNETIC  FIELD  GRADIENT 


The  magnetic  field  gradients  due  to  the  electrode  distributed  currents  can  be  expressed  in 
terms  of  the  r  and  z  derivatives  of  the  magnetic  field  expressions  already  obtained  in  Equations 
(34),  (35),  (77),  (80),  and  (99).  First 


and 


H‘  =  -Hi  sin  0  =  "Hi 


(103) 
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He=H‘cosB=-H‘ 


(104) 


Then,  for  the  diagonal  gradients 


dHex 

dx 


and 


xy  y  dm  _  Xy 
r3  r  dx  r 2 


dHt  dm 


m  dm 


Ig  UIlQ 

r  dr 


(105) 


(106) 


dy  dx  ’ 

the  latter  following  from  V  •  H*  =  0  and  H‘  =  0  (V  •  H  =  0  and  below,  it  is  shown  that  V  •  H®5  =  0). 
For  the  off-diagonal  gradients,  the  non-zero  terms  are 

dm  y  dm 


and 


dz 


r  dz 


dm  xdm 


dm 

dy 


'>  _x _ 

dz  r  dz 


x‘„,  y2dm 


-2  dr 


=h  -m+ 


(107) 


(108) 


(109) 


dmy 

dx 


y2__.  x2dm 


■2  dr 


'-m+- 


(110) 


CONDUCTING  BOTTOM 


Region  1 

From  Equation  (34) 
dHl 


_ ie  _  ^ie  +  I 

dr  r  47C 1 


2nd  +  z+  h 


z+h 


;+  1  Qn{~ 


2nd  -z-h 


'i(r2  +  (z+h?f  "  =  1  yj(r2  +  (2nd  —z  ~h)2f 


2nd  —  z+  h 


2nd  +  z~h 


VO r2  +  {2nd  +  z+hff  V  (r2  +  (2nd  -  z  +  hff  ^ (r2  +  (2nd  +  z  - h)2f 


}]• 


(HI) 
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and 


BHk=J_[ 


( z+hf 


_ "  =  _ [ _ ^  —l - - - j.  y  Qn{ - - — - — - — - - 

&  4nr  yj(r2  +  (z  +  hff  y/r2  +  (z  +h)2  «=»  yj(r2  +  (2nd  -z-hff 

(2nd  +  z+h)2  ,  (2nd -z+hf  ,  (2 nd+z-hf 


(2nd -z-hf 


(112) 


H - —  - . .  =  -■+  p  ■  ' - - - - - =  H - .  ■  - - 

y(r2  +  (2nd  +z+hff  ^(r2 + (2nd -z+hf f  ^(r2  +  (2nd  +  z-hff 

_ 1 _ 1 _ 1 _ 1 

'\r2  + (2nd -z-hf  *v }r2  +  (2nd  +z  +hf  y lr2+(2nd -z  +hf  yjr2  +  (2nd  +z -h) 


fH 


Region  2 

From  Equation  (35) 


and 


a#2e=  m,  I 

dr  -  +^-{ 


z-h 


+(l-Q)lQn\ 

n  =  0 


r  4lt  yj(r2  +  (z-hff 


2nd  +  z-h  2nd  +  z+h 


yl(r2  +  2nd  +z-h2f  (r2 + 2nd  +  z  +  h2f 


(z-hf 


+d-Q)lQn\ 

n  =0 


ZHU  I 


dz  4nr  yj(r2  +  (z  -hff 


(2nd + z-hf  +  (2nd  +  z+hf 


yl(r2  +  2nd+z-h2f  > l(r2+2nd +z+h2f 


(113) 


(114) 


NON-CONDUCTING  BOTTOM  (Q=l) 
Region  1 

From  Equation  (77) 


dffie  /  .  1 

8r  47t  r2 


1  + 


z -h 


2 z 


y!r2  +  (z-hf  d\  V(r2  +  (z-/i)2) 


z -h 


7I2  “ 

H — -  E  n/K0(n7ir/c?)[sin{n7t(z -/7)/cf}  +  sin{n7i(z +/i)/<i}]) 
a  »  =  i 


(115) 


and 
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Bh‘b_  /  ^  1 _ (z -hf  2 

dz  4n  r  y/r2  +  (z-/i)2  V(r2  +  (z  -hff  * 

7C2  ~ 

+r — nKjfanrld) [cos{n n(z  -h)/d}  +  cos{/i %{z  +  h )/d}]). 


Region  2 

From  Equation  (80) 


and 


dHL 


26 


dr 


Hia _ /_'  z-h 

r  4K[^(r2  +  (z~h)2f 


dHlo  _/  f  (z-h)2 
dz  ~  4nr[^f(r2  +  (z-h)2f 


(116) 


(117) 


(118) 


PERFECTLY  CONDUCTING  BOTTOM  (Q=-l,  REGION  1  ONLY) 


From  Equation  (99) 


and 


z -h 


z -h 


271 


dHeie_  I  / 

Br  ~4k(  +  «  + l*Wfl  (119) 


(sin[(/i  +  l/2)7t  (z  -W]  +  sin[(n  +  l/2)7t(z  +  /i  )/</]}) 


a//,', 


16 


3z  4*r  V(r!+(z-/,)2)!  Vr!  +  fe-*)!  rf2-.» 

{cos[(n  +  l/2)7t(z  -/i)/</]  +  cos[(rc  +  1/2)tc(z  +  *)/</]}) 


27tr  “  , 

— TT  ^  (n  +  l/2)/sT,[(n  +  ll2)nr/d]  (120) 


MAGNETIC  FIELD/FIELD  GRADIENT  GENERATED 
BY  A  CABLE  TERMINATED  IN  REGION  1: 
CABLE-ONLY  CONTRIBUTION 


The  magnetic  field/field  gradient  of  a  cable  carrying  a  current/,  with  each  end  shorted  to  the 

medium  in  region  1,  can  be  constructed  as  the  superposition  of  the  magnetic  field/field  gradient  of 

rnr^m  Th  at  ^  Cab- 6  ePd  P0mt  ^  the  Biot-Savart  magnetic  field/field  gradient  of  the  cable 
current.  The  configuration  is  shown  m  Figure  4. 


20 


CSS/TR-94/32 


FIGURE  4.  CURRENT-CARRYING  CABLE/ELECTRODE  CONFIGURATION 


The  electrode  contributions  to  the  magnetic  field  are  constructed  by  means  of  Equations  (34) 
through  (36)  with  appropriate  assignments  of  geometry  and  with  current  signs  as  indicated  in  Figure 
4.  The  electrode  contributions  to  the  magnetic  field  gradient  tensor  are  detailed  in  Equations  (105) 
through  (120).  The  contribution  of  the  cable  is  obtained  by  means  of  a  line  integral  along  the  cable. 
In  practice,  the  line  integral  is  approximated  as  a  sum  of  line  integrals  over  short  linear  segments. 
A  general  representation  is  now  given  of  the  magnetic  field  vector  and  gradient  tensor  for  such  a 
segment. 


BIOT-S AVART  FIELD  AND  GRADIENT  OF  A  STRAIGHT  CABLE  SEGMENT 

The  geometry  for  the  construction  is  shown  in  Figure  5.  Vector  quantities  there  are  indicated 
with  superposed  arrows,  but  in  the  text  they  are  indicated  with  boldface  type.  The  various  vectors 
involved  are 


R„  R2  locate  cable  segment  end  points  (1  — >  2  positive  current  flow) 


R  locates  the  field  point 
R2-Rj  is  along  the  cable  segment 


a 


1  = 


j^2  Rji 
IR2-RJ’ 


Ri  R  R2  R 

u,  =— — — ,  and  u,  =■ 


|Ri-Rf 


IR2-RI 


are  unit  vectors. 


With  these  definitions,  the  general  expression  for  the  magnetic  field  of  the  cable  segment  can  be 
given.  This  expression  is  the  result  of  analytically  performing  the  line  integral 
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„=L  /fl  ,n  n 

47lJ  |  R  —  R'| 3  47c  11 


(121) 


l^-RJ 


=  |;fx(R-R,)  j 


du 


l  (yju2-2l»  (R-Rj)m  +|R-R,|2j3 
over  the  straight  segment,  and  using  some  results  from  vector  geometry.  The  result  is 


1 


R-Ril 


1«  (u2-fl,) 
1  —  (f  •  Uj)2 


and 


H  = 


4tt|R1-Rl 


Uj  xf 


(U2-U!) 

l-(u,.])2 


(122) 


(123) 


z 


FIGURE  5.  CABLE  SEGMENT  GEOMETRY 

The  field  gradient  tensor  is  best  constructed  by  writing  out  Equation  (121 )  in  component  form 
and  taking  the  appropriate  derivatives.  Then 


lRl-R2l 

~  7ZEim«L(^n --^ln)  f  /  r— ; 

471  o  (W 


du 


-  21  •  (R  -  Rj)w  + 1 R  -  Rj| 2]3 


where  the  Levi-Civita  tensor  emn  has  the  properties 


(124) 
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eim« =  1  i,m,n  an  even  permutation  of  1 ,2,3 
£<mn  =  i,m,n  an  odd  permutation  of  1,2,3 
eimn =  0  otherwise. 


This  gives  the  gradient  tensor  components 


dH  J 

3^ = ^  {e, J„  V,  -  3[f  X  (R  -  R,)].  K/tj  -  Rlt)l2  - !/,]} 


where 


IRj-Rjl 


■■  J 


du 


l  (V«2-2f-(R-R1)H+|R-R1|2f 


and 


'■  J 


duu 


(V«2-21*(R-R,)u+|R-R1|2J5’ 


Explicitly,  define 


3= 


1 


l^-RHl-CI-fi,)2] 


2r*(a2-6,) 


l-a-fl,)2 


+1*  (a^-Uj) 


1  a 


12 


|R,-R|2  | R2 - R| 


with 


1  Ri  -  R| 
12  |R2-R| 


Otn  — 


Then 


2  3|  Ri  -  R| 2 

/  = _ ^ _ 

3  3|Rj-R| 

and  the  gradient  tensor  components  can  all  be  written  in  the  compact  form 


(125) 


(126) 

(127) 

(128) 

(129) 

(130) 

(131) 

(132) 

(133) 
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Here,  it  is  noted  that 


dH , 


3^=^  -  a  x  Oi).  (ijx+  uuM 


(134) 


v  •  H  =  -£  {(f  X  Qj)  •  ff?c+  (t x  fij)  •  M  =  o.  (135) 

GRADIENT  MEASUREMENTS:  HOUSING  EFFECTS 


The  presence  of  a  distributed  current  causes  the  undisturbed  gradient  tensor  to  be  asymmetric 
via  the  relation  V  x  H  =  J.  However,  any  practical  gradiometer  will  be  housed  in  a  container  that 
excludes  the  current  from  the  measurement  volume.  This  will  have  a  profound  effect  on  the 
measurement,  in  that  it  may  modify  the  magnetic  field  in  the  enclosure  relative  to  the  original  field, 
and  it  forces  the  measured  gradient  tensor  to  be  symmetric,  possibly  in  a  way  that  may  be  difficult 
to  determine. 

Joseph,  et.  al.,u’3  have  analyzed  the  effects  of  enclosures  to  determine  the  relationship  between 
the  measured  field  vector  and  gradient  tensor  components  and  those  that  exist  in  the  absence  of  the 
enclosure.  They  have  given  detailed  results  for  spherical  and  axisymmetric  enclosures  under  the 
assumption  that  there  are  no  other  boundaries,  and  that  the  undisturbed  current  density  J°  is  uniform 
over  the  volume  occupied  by  the  enclosure. 

SPHERICAL  ENCLOSURE 

For  a  spherical  enclosure,  the  interior  field  has  the  form 

H  =  H°  +  H'  (136) 

where  H°  is  the  field  in  the  absence  of  the  enclosure,  and  H'  is  the  distortion  field,  which  is  related 
to  the  current  J°  via 


H'  =  ^(rxJ°)  (137) 

From  this  expression,  it  is  clear  that  the  enclosure  has  no  effect  on  the  field  at  the  center  of 
the  sphere.  Away  from  the  center,  there  is  a  nonzero  correction  term.  The  expression  also  can  be 
used  to  show  that  the  effect  of  the  enclosure  on  the  gradient  tensor  is  to  select  the  symmetric  part 
of  the  tensor;  that  is. 


(Gj  +  Gj) 


(138) 


an  expression  that  holds  everywhere  within  the  enclosure.  Thus,  the  measured  gradient  tensor  may 
be  easily  calculated  in  terms  of  the  undistorted  tensor. 
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AXISYMMETRIC  ENCLOSURE 

In  the  axisymmetric  case,  the  enclosure  affects  the  field  in  the  following  way.  On  the  symmetry 
axis,  the  axial  distortion  field  is  zero,  while  the  field  perpendicular  to  (athwart)  the  axis  has  nontrivial 
corrections.  Off  the  symmetry  axis,  both  components  of  the  field  have  nontrivial  corrections. 

For  the  gradient  tensor,  on  the  symmetry  axis,  the  diagonal  components  are  undistorted,  and 
those  involving  derivatives  with  respect  to  the  coordinates  perpendicular  to  the  symmetry  axis  are 
symmetrized  as  in  the  spherical  case.  Those  components  involving  derivatives  with  respect  to  the 
axial  coordinate  have  corrections  that  depend  on  both  the  ambient  current  density  and  the  axial 
enclosure  profile,  and  these  corrections  may  be  as  large  as  the  undisturbed  gradients. 


GENERAL  ENCLOSURE 

In  the  general  nonsymmetric  enclosure  case,  all  the  gradient  tensor  components  have  nontrivial 
corrections  depending  on  the  current  density  vector  and  die  enclosure  shape.  The  determination  of 
the  corrections  involves  the  solution  of  the  system  of  equations3 


hW/svVx(j°-oVY) 

(139) 

vy=o 

(140) 

f  0 

oo 

(141) 

aft  •  V<f>'  =  ft  •  J°|  s 

(142) 

where  <J>  is  the  potential  due  to  the  presence  of  the  enclosure  in  the  ambient  (uniform)  current 
distribution,  a  is  the  medium  conductivity,  and  ft  is  the  outward  normal  to  the  surface  of  the 
enclosure. 


The  computation  of  enclosure  corrections  is  beyond  the  scope  of  this  report.  It  will  be  limited 
to  the  case  of  a  spherical  enclosure.  For  this  case  the  measured  gradient  tensor  from  the 
electrode-cable  combination  will  have  the  form 


Gs  = 


|  (g;+g;+g;+g;) 
1 


-(g;+g;+g;+g;)  -(g;+g;+g;) 


(g;-g;> 


i 


2(g;+g;+g;)  ^(g;+g;+g;) 


-(g;+g;+g;) 

-(G*+Gp 


(143) 
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APPENDIX  A 

TRANSERA  HIGH  TECH  BASIC  CODES  FOR  REGION  1 


A-l/A-2 
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CURRENT  DENSITY  CODE 


10  SUB  Jelectrode(R(*),Re(*),Db,S,Sb,I,Jay(*),Err) 

20  X=R(1) 

30  Y=R(2) 

40  Z=R(3) 

50  Xe=Re(l) 

60  Ye=Re(2) 

70  Ze=Re(3) 

80  Dx=X-Xe 
90  Dy=Y-Ye 

100  IF  SQR((Dx*Dx+Dy  *Dy)/(X*X+Y*  Y +Xe  *Xe+Y  e  *  Ye))<  1  .E- 1 2  THEN 

110  Jay(l)=0 

120  Jay(2)=0 

130  Jay(3)=0 

140  GOTO  1560 

150  END  IF 

160  ! 

170  !This  program  calculates  the  current  density  within  the  upper 
180  llayer  of  a  two-layer  conductor,  bounded  above  by  air,  due  to 

190  !an  electrode  injecting  a  current  I  into  that  layer.  X,Y,  and 
200  !Z  are  the  coordinates  (meters)  of  the  field  point  in  a  coordi- 
210  !nate  system  whose  origin  is  at  the  surface.  Xe,  Ye,  and  Ze  are 
220  !the  coordinates  (meters)  of  the  electrode.  Db  is  the  depth  to 
230  !the  conductor-conductor  boundary  below  the  air-conductor  bound- 

240  !ary.  S  is  the  conductivity  (Siemens/meter)  of  the  upper  conduct- 
250  ling  layer  and  Sb  is  that  for  the  lower  layer(for  Sb=0,  corre- 
260  Isponding  to  an  insulating  bottom,  or  Sb<0,  a  perfectly  conducting 
270  ! bottom,  alternate  series  expressions  are  used).  Jay(*)  is  the 

280  !  three-component  field  vector  in  amperes/meter  A2.  Err  is  the  size 

290  !of  the  summand  relative  to  the  sum  when  summation  ceases.  Note 
300  Ithat  I  is  positive  if  current  enters  the  medium  from  the  electrode 

310  land  negative  otherwise. 

320  ! 

330  OPTION  BASE  1 
340  IF  Sb>0  THEN 
350  Q=(S-Sb)/(S+Sb) 

360  ELSE 

370  IF  Sb=0  THEN 

380  Q=1 

390  ELSE 

400  Q=-l 

410  END  IF 

420  END  IF 

430  X2=(X-Xe)  *(X-Xe) 

440  Y2=(Y -Ye)*(Y-Ye) 

450  R0=SQR(X2+Y  2) 

460  R03=R0*R0*R0 
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470  Rm=SQR(R0*R0+(Z-Ze)*(Z-Ze)) 

480  Rm3=Rm*Rm*Rm 

490  Rp=S  QR(R0  *R0+(Z+Ze)  *  (Z+Ze)) 

500  Rp3=Rp*Rp*Rp 

510  IF  ABS(Q)<1  THEN 

520  Smxy=l/Rm3+1/Rp3 

530  Smz=(Z-Ze)/Rm3+(Z+Ze)/Rp3 

540  Qn=Q 

550  N=1 

560  Zmm=2*N*Db-Z-Ze 
570  Zpp=2*N*Db+Z+Ze 
580  Zmp=2*N*Db-Z+Ze 
590  Zpm=2*N*Db+Z-Ze 
600  Rmm=SQR(R0*R0+Zmm*Zmm) 

610  Rpp=SQR(R0*R0+Zpp*Zpp) 

620  Rmp=SQR(R0*R0+Zmp*Zmp) 

630  Rpm=SQR(R0*R0+Zpm*Zpm) 

640  Rmm3=Rmm*Rmm*Rmm 
650  Rpp3=Rpp*Rpp*Rpp 
660  Rmp3=Rmp*Rmp*Rmp 
670  Rpm3=Rpm*Rpm*Rpm 

680  Intz=Qn*(-Zmm/Rmm3+Zpp/Rpp3-Zmp/Rmp3+Zpm/Rpm3) 
690  Intxy=Qn*(l/Rmm3+l/Rpp3+l/Rmp3+l/Rpm3) 

700  Smz=Smz+Intz 

710  Smxy=Smxy+Intxy 

720  IF  N>4  THEN 

730  IF  ABS(Intz/Smz)<Err  THEN 

740  IF  ABS(Intxy/Smxy)<Err  THEN 

750  GOTO  820 

760  END  IF 

770  END  IF 

780  END  IF 

790  N=N+1 

800  Qn=rQ*Qn 

810  GOTO  560 

820  Jay(3)=I*Smz/4/PI 

830  Jay(l)=(X-Xe)*I*Smxy/4/PI 

840  Jay(2)=(Y-Ye)*I*Smxy/4/PI 

850  ELSE 

860  Zm=Z-Ze 

870  Zp=Z+Ze 

880  Smxy=0 

890  Smz=0 

900  Fxy=0 

910  Fz=0 

920  IF  Q=1  THEN 

930  N=1 

940  CALL  K0(N*PI*R0/Db,K0) 

950  CALL  Kl(N*PI*R0/Db,Kl) 

960  Intxy=N*K  1  *(COS(N*PI*Zm/Db)+COS(N*PI*Zp/Db)) 

970  Intz=N*KO*(SIN(N*PI*Zm/Db)+SIN(N*PI*Zp/Db)) 

980  Smxy=Smxy+Intxy 
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990  Smz=Smz+Intz 

1000  IF  N>4  THEN 

1010  IF  ABS(Smxy)=0  THEN 

1020  Fxy=l 

1030  ELSE 

1040  IF  ABS(Intxy)/ABS(Smxy)<Err  THEN 

1050  Fxy=l 

1060  END  IF 

1070  END  IF 

1080  IF  ABS(Smz)=0  THEN 

1090  Fz=l 

1100  ELSE 

1110  IF  ABS(Intz)/ABS(Smz)<Err  THEN 

1120  Fz=l 

1130  END  IF 

1140  END  IF 

1 150  IF  Fxy+Fz>l  THEN  GOTO  1 190 

1160  N=N+1 

1170  GOTO  940 

1180  END  IF 

1190  Jay(3)=I/2/Db/Db*Smz 

1200  Jfac=I/(2*PI*R0*R0*Db)*(l+PI*R0/Db*Smxy) 

1210  Jay(  1  )=(X-Xe)*  Jfac 

1220  Jay(2)=(  Y  -  Y  e)*  Jfac 

1230  ELSE 

1240  N=0 

1250  CALL  K0((N+.5)*PI*R0/Db,K0) 

1260  CALL  Kl((N+.5)*PI*R0/Db,Kl) 

1270  Intxy=(N+.5)*Kl  *(COS((N+.5)*PI*Zm/Db)+COS((N+.5)*PI*Zp/Db)) 

1280  Intz=(N+.5)*K0*(SIN((N+.5)*PI*Zm/Db)+SIN((N+.5)*PI*Zp/Db)) 

1290  Smxy=Smxy+Intxy 

1300  Smz=Smz+Intz 

1310  IF  N>4  THEN 

1320  IF  ABS(Smxy)=0  THEN 

1330  Fxy=l 

1340  ELSE 

1350  IF  ABS(Intxy)/ABS(Smxy)<Err  THEN 

1360  Fxy=l 

1370  END  IF 

1380  END  IF 

1390  IF  ABS(Smz)=0  THEN 

1400  Fz=l 

1410  ELSE 

1420  IF  ABS(Intz)/ABS(Smz)<Err  THEN 

1430  Fz=l 

1440  END  IF 

1450  END  IF 

1460  IF  Fxy+Fz>l  THEN  GOTO  1500 

1470  N=N+1 

1480  GOTO  1250 

1490  END  IF 

1500  Jay(3)=I/2/Db/Db*Smz 
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1510  Jfac=l/(2*R0*Db*Db) 
1520  Jay(  1  )=(X-Xe)*  Jfac*Smxy 

1530  Jay(2)=(Y-Ye)*Jfac*Smxy 
1540  END  IF 
1550  END  IF 
1560  SUBEND 


ELECTRODE  MAGNETIC  FIELD  CODE 


10  SUB  HeIectrode(R(*),Re(*),Db,S,Sb,I,Aitch(*),EiT) 

20  X=R(1) 

30  Y=R(2) 

40  Z=R(3) 

50  Xe=Re(l) 

60  Ye=Re(2) 

70  Ze=Re(3) 

80  Dx=X-Xe 
90  Dy=Y-Ye 

100  IF  SQR((Dx*Dx+Dy*Dy)/(X*X+Y*  Y +Xe*Xe+Ye*Ye))<l  .E- 1 2  THEN 
110  Aitch(l)=0 
120  Aitch(2)=0 
130  Aitch(3)=0 
140  GOTO  1080 
150  END  IF 
160  ! 

170  !This  program  calculates  the  magnetic  field  within  the  upper 
180  Mayer  of  a  two-layer  conductor,  bounded  above  by  air,  due  to 
190  !an  electrode  injecting  a  current  I  into  that  layer.  X,Y,  and 
200  !Z  are  the  coordinates  (meters)  of  the  field  point  in  a  coordi- 
210  Inate  system  whose  origin  is  at  the  surface.  Xe,  Ye,  and  Ze  are 
220  !the  coordinates  (meters)  of  the  electrode.  Db  is  the  depth  to 
230  !the  conductor-conductor  boundary  below  the  air-conductor  bound- 

240  !ary.  S  is  the  conductivity  (Siemens/meter)  of  the  upper  conduct- 
250  !ing  layer  and  Sb  is  that  for  the  lower  layer(for  Sb=0,  corre- 
260  Isponding  to  an  insulating  bottom,  or  Sb<0,  a  perfectly  conducting 
270  Ibottom,  alternate  series  expressions  are  used).  Aitch(*)  is 

280  !the  three-component  field  vector(ampe res/meter).  Err  is  the  size 

290  lof  the  summand  relative  to  the  sum  when  summation  ceases.  Note 
300  !that  I  is  positive  if  current  enters  the  medium  from  the  electrode 

310  land  negative  otherwise. 

320  ! 

330  OPTION  BASE  1 

340  IF  Sb>0  THEN 
350  Q=(S-Sb)/(S+Sb) 

360  ELSE 

370  IF  Sb=0  THEN 

380  Q=1 

390  ELSE 

400  Q=-l 

410  END  IF 
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420  END  IF 

430  X2=(X-Xe)*(X-Xe) 

440  Y2=(Y -Ye)*(Y-Ye) 

450  R0=SQR(X2+Y2) 

460  R03=R0*R0*R0 

470  Rp=SQR(R0*R0+(Z+Ze)*(Z+Ze)) 

480  R3=Rp*Rp*Rp 
490  IF  ABS(Q)<1  THEN 
500  Sm=l-(Z+Ze)/Rp 
510  Qn=Q 
520  N=1 

530  Zmm=2*N*Db-Z-Ze 
540  Zpp=2*N*Db+Z+Ze 
550  Zmp=2*N*Db-Z+Ze 
560  Zpm=2*N*Db+Z-Ze 
570  Rmm=SQR(R0*R0+Zmm*Zmm) 

580  Rpp=SQR(R0*R0+Zpp*Zpp) 

590  Rmp=SQR(R0*R0+Zmp*Zmp) 

600  Rpm=SQR(R0*R0+Zpm*Zpm) 

610  Int=Qn*(Zmm/Rmm-Zpp/Rpp+Zmp/Rmp-Zpm/Rpm) 

620  Sm=Sm+Int 

630  IF  N>2  THEN 

640  IF  ABS(Int/Sm)<Err  THEN 

650  Sm=I*Sm/(4*PI*R0) 

660  GOTO  1050 

670  END  IF 
680  END  IF 
690  N=N+1 

700  Qn=Q*Qn 
710  GOTO  530 
720  ELSE 
730  Sm=0 
740  Zm=Z-Ze 
750  Zp— Z+Ze 

760  Rm=SQR(R0*R0+Zm*Zm) 

770  IF  Q=1  THEN 
780  N=1 

790  CALL  Kl(N*PI*R0/Db,Kl) 

800  Int=Kl*(SIN(N*PI*Zm/Db)+SIN(N*PI*Zp/Db)) 

810  Sm=Sm+Int 

820  IF  N>=2  THEN 

830  IF  ABS(Int)/ABS(Sm)<Err  THEN  880 

840  ELSE 

850  N=N+1 

860  GOTO  790 

870  END  IF 

880  Sm=l-2*Z/Db+Zm/Rm+PI*R0/Db*Sm 

890  Sm=I*Sm/(4*PI*R0) 

900  ELSE 
910  N=0 

920  CALL  Kl((N+.5)*PI*R0/Db,Kl) 

930  Int=Kl  *(SIN((N+.5)*PPZm/Db)+SIN((N+.5)*PI*Zp/Db)) 
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940  Sm=Sm+Int 

950  IF  N>=2  THEN 

960  IF  ABS(Int)/ABS(Sra)<Err  THEN  1010 

970  ELSE 

980  N=N+1 

990  GOTO  920 

1000  END  IF 

1010  Sm=l+Zm/Rm-2*R0/Db*Sm 
1020  Sm=I*Sm/(4*PI*R0) 

1030  END  IF 
1040  END  IF 

1050  Aitch(  1  )=(Y  e- Y)  *Sm/R0 
1060  Aitch(2)=(X-Xe)  *Sm/R0 
1070  Aitch(3)=0 
1080  SUBEND 


ELECTRODE  MAGNETIC  FIELD  GRADIENT  CODE 


10 

20 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 


SUB  GeIectrode(R(*),Re(*),Db,S,Sb,I,G(*),Err) 

OPTION  BASE  1 
DIM  H(3) 

X=R(1) 

Y=R(2) 

Z=R(3) 

Xe=Re(l) 

Ye=Re(2) 

Ze=Re(3) 

t 

• 

!This  program  calculates  the  magnetic  field  gradient  within  the  upper 
Mayer  of  a  two-layer  conductor,  bounded  above  by  air,  due  to 
Jan  electrode  injecting  a  current  I  into  that  layer.  X,Y,  and 
!Z  are  the  coordinates  (meters)  of  the  field  point  in  a  coordi- 
Jnate  system  whose  origin  is  at  the  surface.  Xe,  Ye,  and  Ze  are 
•the  coordinates  (meters)  of  the  electrode.  Db  is  the  depth  to 
Jthe  conductor-conductor  boundary  below  the  air-conductor  bound- 
.ary.  S  is  the  conductivity  (Siemens/meter)  of  the  upper  conduct- 
ling  layer  and  Sb  is  that  for  the  lower  Iayer(for  Sb=0,  corre- 
Jsponding  to  an  insulating  bottom,  or  Sb<0,  a  perfectly  conducting 
Jbottom,  alternate  series  expressions  are  used).  G(*)  is  the 
!3x3  magnetic  gradient  tensor(amperes/meterA2).  Err  is  the  size 
!°f  the  summand  relative  to  the  sum  when  summation  ceases.  Note 
•that  I  is  positive  if  current  enters  the  medium  from  the  electrode 
land  negative  otherwise. 

t 

IF  Sb>0  THEN 
Q=(S-Sb)/(S+Sb) 

ELSE 

IF  Sb=0  THEN 

Q=l 

ELSE 
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330  Q=-l 

340  END  IF 

350  END  IF 

360  X2=(X-Xe)*(X-Xe) 

370  Y2=(Y -Ye)*(Y-Ye) 

380  Xy=(X-Xe)*(Y  -Ye) 

390  R0=SQR(X2+Y  2) 

400  R02=R0*R0 

410  R03=R0*R0*R0 

420  Rp=SQR(R0*R0+(Z+Ze)  *(Z+Ze)) 

430  R3=Rp*Rp*Rp 

440  IF  ABS(Q)<1  THEN 

450  Smtr=(Z+Ze)/R3 

460  Smtz=(Z+Ze)*(Z+Ze)/R3-  1/Rp 

470  Qn=Q 

480  N=1 

490  Zmm=2*N*Db-Z-Ze 

500  Zmm2=Zmm*Zmm 

510  Zpp=2*N*Db+Z+Ze 

520  Zpp2=Zpp*Zpp 

530  Zmp=2*N*Db-Z+Ze 

540  Zmp2=Zmp*Zmp 

550  Zpm=2*N*Db+Z-Ze 

560  Zpm2=Zpm*Zpm 

570  Rmm=SQR(R0*R0+Zmm*Zmm) 

580  Rmm3=Rmm*Rmm*Rmm 
590  Rpp=SQR(R0*R0+Zpp*Zpp) 

600  Rpp3=Rpp*Rpp*Rpp 

610  Rmp=SQR(R0*R0+Zmp*Zmp) 

620  Rmp3=Rmp*Rmp*Rmp 
630  Rpm=SQR(R0*R0+Zpm*Zpm) 

640  Rpm3=Rpm*Rpm*Rpm 

650  Inttr=Qn*(-Zmm/Rmm3+Zpp/Rpp3-Zmp/Rmp3+Zpm/Rpm3) 

660  Intt2=Qn*(Zmm2/Rmm3+Zpp2/Rpp3+Zmp2/Rmp3+Zpm2/Rpra3) 
670  Inttz=Inttz+Qn*(-  1/Rmm-  1/Rpp-  1/Rmp-  1/Rpm) 

680  Smtr=Smtr+Inttr 

690  Smtz=Smtz+Inttz 

700  IF  N>2  THEN 

710  IF  ABS(Inttr/Smtr)<Err  THEN 

720  IF  ABS(Inttz/Smtz)<Err  THEN 

730  Smtr=I*Smtr/(4*PI) 

740  Smtz=I*Smtz/(4*PI*R0) 

750  GOTO  1290 

760  END  IF 

770  END  IF 

780  END  IF 
790  N=N+1 

800  Qn=Q*Qn 
810  GOTO  490 
820  ELSE 
830  Smtr=0 
840  Smtz=0 
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850  Zm=Z-Ze 
860  Zp=Z+Ze 

870  Rm=SQR(R0*R0+Zm*Zm) 

880  Rm3=Rm*Rm*Rm 
890  IF  Q=1  THEN 
900  N=1 

910  CALL  K0(N*PI*R0/Db,K0) 

920  Inttr=N*KO*(SIN(N*PI*Zm/Db)+SIN(N*PI*Zp/Db)) 

930  Sratr=Smtr+Inttr 

940  CALL  Kl(N*PI*R0/Db,Kl) 

950  Inttz=N*Kl*(COS(N*PI*Zm/Db)+COS(N*PI*Zp/Db)) 

960  Smtz=Smtz+Inttz 

970  IF  N>=2  THEN 

980  IF  ABS(Inttr)/ABS(Smtr)<Err  THEN 

990  IF  ABS(Inttz)/ABS(Smtz)<Err  THEN 

1000  GOTO  1060 

1010  END  IF 

1020  END  IF 

1030  END  IF 

1040  N=N+1 

1050  GOTO  910 

1060  Smtr=I*(-(  l+Zm/Rm-2*Z/Db)/R02-Zm/Rm3+PI*PI/Db/Db*Smtr)/4/PI 

1070  Smtz=I*(-2/Db-Zm*Zm/Rm3+l/Rm+PI*PI*R0*Smtz/Db/Db)/(4*PI*R0) 
1080  ELSE 
1090  N=0 

1 100  CALL  K0((N+.5)*PI*R0/Db,K0) 

1 1 10  Inttr=(N+.5)*K0*(SIN((N+.5)*PI*Zm/Db)+SIN((N+.5)*PI*Zp/Db)) 

1120  Smtr=Smtr+Inttr 

1130  CALL  Kl((N+.5)*PI*R0/Db,Kl) 

1 140  Inttz=(N+.5)*Kl*(COS((N+.5)*PI*Zm/Db)+COS((N+.5)*PI*Zp/Db)) 

1150  Smtz=Smtz+Inttz 

1160  IF  N>=2  THEN 

1170  IF  ABS(Inttr)/ABS(Smtr)<Err  THEN 

1180  'IF  ABS(Inttz)/ABS(Smtz)<Err  THEN 

1190  GOTO  1250 

1200  END  IF 

1210  END  IF 

1220  END  IF 

1230  N=N+1 

1240  GOTO  1100 

1250  Smtr=-I*(Zm/Rm/R02+Zm/Rm3+2*PI*Srntr/Db/Db)/4/PI 
1260  Smtz=I*(-Zm*Zm/Rm3+l/Rm-2*PI*R0*Smtz/Db/Db)/(4*PI*R0) 

1270  END  IF 
1280  END  IF 

1290  CALL  Helectrode(R(*),Re(*),Db,S,Sb,I,H(*),Err) 

1 300  Ht=(X-Xe)*H(2)/R0-(Y -Ye)*H(  1  )/R0 
1310  IF  ABS(Q)<1  THEN 
1320  Smtr=Smtr-Ht/R0 
1330  END  IF 

1340  G(  1 , 1  )=Xy  *(Ht-R0*Smtr)/R03 
1350  G(2,2)=-G(l,l) 

1360  G(3,3)=0 
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1370  G(3,l)=0 
1380  G(3,2)=0 

1390  G(l,2)=-(X2*Ht+R0*Y2*Smtr)/R03 
1400  G(2, 1  )=(Y2*Ht+R0*X2*Smtr)/R03 
1410  G(1,3)=-(Y -Ye)*Smtz/R0 
1420  G(2,3)=(X-Xe)*Smtz/R0 
1430  SUBEND 


BIOT-SA  VART  MAGNETIC  FIELD  CODE 


10  SUB  Hbiotsavrt(I,Rl(*),R2(*),R(*),H(*)) 

20  OPTION  BASE  1 

30  DIM  X 1  (3),X2(3),X  12(3),L(3),U  1  (3),U2(3) 

40  ! 

50  !This  program  calculates  the  magnetic  field  H(*)  (amperes/meter) 
60  !at  the  point  R(*)  (meters),  due  to  a  straight  wire  segment 
70  !whose  end  points  are  Rl(*)  and  R2(*)  (meters),  and  in  which 
80  !the  current  flows  from  Rl(*)  to  R2(*). 

90  ! 

100  X1(1)=R1(1)-R(1) 

110  X 1  (2)=R  1  (2)-R(2) 

120  X 1  (3)=R1  (3)-R(3) 

130  X2(1)=R2(1)-R(1) 

140  X2(2)=R2(2)-R(2) 

150  X2(3)=R2(3)-R(3) 

160  X12(1)=R2(1)-R1(1) 

170  X 1 2(2)=R2(2)-R  1  (2) 

180  X 1 2(3)=R2(3)-R  1(3) 

190  Rlm=SQR(DOT(Xl,Xl)) 

200  R2m=SQR(DOT(X2,X2)) 

210  R12m=SQR(DOT(X12,X12)) 

220  MAT  L=  (l/R12m)*X12 
230  MAT  Ul=  (l/Rlm)*Xl 
240  MAT  U2=  ( l/R2m)*X2 
250  N=DOT(L,U  1  )-DOT(L,U2) 

260  D=l-DOT(L,Ul)*DOT(L,Ul) 

270  Fac=(I*N)/(4*PI*Rlm*D) 

280  H(  1  )=L(2)  *U  1  (3)-L(3)  *U  1  (2) 

290  H(2)=L(3)  *U  1  ( 1  )-L(  1 )  *U  1  (3) 

300  H(3)=L(  1 )  *U  1  (2)-L(2)  *U  1  ( 1 ) 

310  MAT  H=  (Fac)*H 
320  SUBEND 


BIOT-SA  YART  MAGNETIC  FIELD  GRADIENT  CODE 


10  SUB  Gbiotsavrt(I,Rl(*),R2(*),R(*),G(*)) 
20  OPTION  BASE  1 

30  DIM  X 1  (3),X2(3)^X  1 2(3),L(3),U  1  (3),U2(3) 
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40  ! 

50  !This  program  calculates  the  magnetic  field  gradient  tensor 
60  !G(*)  (ampere/meterA2)  at  the  point  R(*)  (meters),  due  to  a 

70  Istraight  wire  segment  whose  end  points  are  Rl(*)  and  R2(*) 

80  !(meters)  and  in  which  positive  current  I  flows  from  Rl(*) 

90  !to  R2(*). 

100  ! 

110  X1(1)=R1(1)-R(1) 

120  X1(2)=R1(2)-R(2) 

130  X1(3)=R1(3)-R(3) 

140  X2(1)=R2(1)-R(1) 

150  X2(2)=R2(2)-R(2) 

160  X2(3)=R2(3)-R(3) 

170  X12(1)=R2(1)-R1(1) 

180  X  12(2)=R2(2)-R  1  (2) 

190  X 1 2(3)=R2(3  )-R  1  (3) 

200  R 1  m=SQR(DOT  (X 1  ,X  1 )) 

210  R2m=SQR(DOT(X2,X2)) 

220  R 1 2m=SQR(DOT  (X 1 2,X  1 2) ) 

230  MAT  L=(l/R12m)*X12 
240  MAT  Ul=(l/Rlm)*Xl 
250  MAT  U2=(l/R2m)*X2 
260  N=DOT(L,Ul)-DOT(L,U2) 

270  D=l-DOT(L,Ul)*DOT(L,Ul) 

280  1 1  =-N/D/R  1  m/R  1  m 

290  A12=Rlm/R2m 

300  I2=(-2*N/D+A  1 2*  A 1 2*DOT(L,U2)-DOT(L,U  1  ))/D/R  1  m/R  1  m 

310  I3=l/R  lm/R  1  m-A  1 2/R2m/R2m-DOT(L,U  1  )*I2 

320  Gfac=I/4/PI 

330  G(  1 , 1  )=Gfac*(-L(2)*U  1  (3)+L(3)*U  1  (2))*(L(  1  )*I3+U  1  ( 1  )*I2) 

340  G(2,2)=Gfac*(-L(3)*U  1  ( 1  )+L(  1  )*U  1  (3))*(L(2)*I3+U  1  (2)*I2) 

350  G(3,3)=-G(  1 ,1  )-G(2,2) 

360  G(1 ,2)=Gfac*(-L(3)*I  1  +(-L(2)*U  1  (3)+L(3)*U  1  (2))*(L(2)*I3+U  1  (2)*I2» 

370  G(2,l  )=Gfac*(L(3)*I  1  +(-L(3)*U  1  ( 1  )+L(  1  )*U  1  (3))*(L(  1  )*I3+U  1  ( 1  )*I2» 

380  G(l,3)=Gfac*(L(2)*Il+(-L(2)*Ul(3)+L(3)*Ul(2))*(L(3)*I3+Ul(3)*I2)) 
390  G(3,l)=Gfac*(-L(2)*Il+(-L(l)*Ul(2)+L(2)*Ul(l))*(L(l)*I3+Ul(l)*I2» 
400  G(2,3)=Gfac*(-L(  1  )*I1  +(-L(3)  *U  1  ( 1  )+L(  1  )*U  1  (3))*(L(3)*I3+U  1  (3)*I2)) 

4 1 0  G(3,2)=Gfac*(L(  1  )*I  1  +(-L(  1  )*U1  (2)+L(2)*U  1  ( 1  ))*(L(2)*I3+U  1  (2)*I2)) 


MODIFIED  BESSEL  FUNCTION  CODES 


10  SUB  I0(X,I0) 

20  ! 

30  Icalculates  the  modified  bessel  function 
40  ! 

50  OPTION  BASE  1 

60  DIM  A(9),T(9) 

70  MAT  A=(0) 

80  MAT  T=(0) 


I0(x)  for  real  x>=0 
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90  U=X/3.75 

100  IF  U<=1  THEN 
110  A(l)=l 
120  A(2)=3.5 156229 

130  A(3)=3.0899424 

140  A(4)=1.2067492 

150  A(5)=.2659732 

160  A(6)=.0360768 

170  A(7)=.0045813 

180  T(l)=l 
190  T(2)=U*U 

200  T(3)=T(2)*T(2) 

210  T(4)=T(2)*T(3) 

220  T(5)=T(3)*T(3) 

230  T(6)=T(2)*T(5) 

240  T(7)=T(4)*T(4) 

250  IO=DOT(A,T) 

260  ELSE 
270  A(l)=.39894228 

280  A(2)=.0 1328592 

290  A(3)=. 002253 19 

300  A(4)=-.00157565 

310  A(5)=.009 16281 

320  A(6)=-.02057706 

330  A(7)=.02635537 

340  A(8)=-.01647633 

350  A(9)=.00392377 

360  T(l)=l 
370  T(2)=l/U 

380  T(3)=T(2)*T(2) 

390  T(4)=T(2)*T(3) 

400  T(5)=T(3)*T(3) 

410  T(6)=T(2)*T(5) 

420  T(7)=T(4)*T(4) 

430  T(8)=T(2)*T(7) 

440  T(9)=T(5)*T(5) 

450  IO=EXP(X)/SQR(X)*DOT(A,T) 

460  END  IF 
470  SUBEND 

480  SUB  I1(X,I1) 

490  ! 

500  Icalculates  the  modified  bessel  function  Il(x)  for  real  x>=0 
510  ! 

520  OPTION  BASE  1 
530  DIM  A(9),T(9) 

540  MAT  A=(0) 

550  MAT  T=(0) 

560  U=X/3.75 

570  IF  U<=1  THEN 
580  A(l)=.5 

590  A(2)=. 87890594 
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600  A(3)=.5 1498869 

610  A(4)=.  15084934 

620  A(5)=.02658733 

630  A(6)=.00301532 

640  A(7)=.00032411 

650  T(l)=l 
660  T(2)=U*U 

670  T(3)=T(2)*T(2) 

680  T(4)=T(2)*T(3) 

690  T(5)=T(3)*T(3) 

700  T(6)=T(2)*T(5) 

710  T(7)=T(4)*T(4) 

720  Il=X*DOT(A,T) 

730  ELSE 
740  A(1  >=.39894228 

750  A(2)=-.03988024 

760  A(3)=-.00362018 

770  A(4)=.00163801 

780  A(5)=-.01031555 

790  A(6)=.02282967 

800  A(7)=-.02895312 

810  A(8)=.0 1787654 

820  A(9)=-.00420059 

830  T(l)=l 
840  T(2)=l/U 

850  T(3)=T(2)*T(2) 

860  T(4)=T(2)*T(3) 

870  T(5)=T(3)*T(3) 

880  T(6)=T(2)*T(5) 

890  T(7)=T(4)*T(4) 

900  T(8)=T(2)*T(7) 

910  T(9)=T(5)*T(5) 

920  Il=EXP(X)/SQR(X)*DOT(A,T) 

930  END  IF 
940  SUBEND 
950  SUB  K0(X,K0) 

960  ! 

970  ! calculates  the  modified  bessel  function  K0(x)  for  real  x>=0 

980  ! 

990  OPTION  BASE  1 
1000  DIM  A(7),T(7) 

1010  MAT  A=(0) 

1020  MAT  T=(0) 

1030  U=X/2 
1040  IF  U<=1  THEN 
1050  A(l)=-.57721566 

1060  A(2)=.42278420 

1070  A(3)=.23069756 

1080  A(4)=.03488590 

1090  A(5)=.00262698 

1100  A(6)=.00010750 

1110  A(7)=.00000740 
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1120  T(l)=l 
1130  T(2)=U*U 

1140  T(3)=T(2)*T(2) 

1150  T(4)=T(2)*T(3) 

1160  T(5)=T(3)*T(3) 

1170  T(6)=T(2)*T(5) 

1180  T(7)=T(4)*T(4) 

1190  CALL  I0(X,I) 

1200  KO=-LOG(U)  *I+DOT  ( A,T) 

1210  ELSE 

1220  A(l)= 1.2533 14 14 

1230  A(2)=-.07832358 

1240  A(3)=.02 189568 

1250  A(4)=-.0 1062446 

1260  A(5)=.00587872 

1270  A(6)=-.00251540 

1280  A(7)=.00053208 

1290  T(l)=l 
1300  T(2)=l/U 

1310  T(3)=T(2)*T(2) 

1320  T(4)=T(2)*T(3) 

1330  T(5)=T(3)*T(3) 

1340  T(6)=T(2)*T(5) 

1350  T(7)=T(4)*T(4) 

1360  KO=EXP(-X)/SQR(X)*DOT(A,T) 

1370  END  IF 
1380  SUBEND 
1390  SUB  K1(X,K1) 

1400  ! 

1410  ! calculates  the  modified  bessel  function  Kl(x)  for  real  x> 
1420  ! 

1430  OPTION  BASE  1 
1440  DIM  A(7),T(7) 

1450  MAT  A=(0) 

1460  MAT  T=(0) 

1470  U=X/2 
1480  IF  U<=1  THEN 
1490  A(l)=l 
1500  A(2)=.  15443 144 

1510  A(3)=-.67278579 

1520  A(4)=-.  18 156897 

1530  A(5)=-.01919402 

1540  A(6)=-.001 10404 

1550  A(7)=-.00004686 

1560  T(l)=l 
1570  T(2)=U*U 

1580  T(3)=T(2)*T(2) 

1590  T(4)=T(2)*T(3) 

1600  T(5)=T(3)*T(3) 

1610  T(6)=T(2)*T(5) 

1620  T(7)=T(4)*T(4) 

1630  CALL  I1(X,I) 
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1640  K 1  =LOG(U)  *I+DOT(  A,T)/X 

1650  ELSE 

1660  A(l)=1.25331414 

1670  A(2)=.23498619 

1680  A(3)=-.03655620 

1690  A(4)=.0 1504268 

1700  A(5)=-.00780353 

1710  A(6)=.00325614 

1720  A(7)=-.00068245 

1730  T(l)=l 
1740  T(2)=l/U 

1750  T(3)=T(2)*T(2) 

1760  T(4)=T(2)*T(3) 

1770  T(5)=T(3)*T(3) 

1780  T(6)=T(2)*T(5) 

1790  T(7)=T(4)*T(4) 

1800  K 1  =EXP(-X)/SQR(X)*DOT(  A,T) 
1810  END  IF 
1820  SUBEND 
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